Nano-indentation of circular graphene flakes 
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Nano-indentation of circular graphene flakes are studied by using Molecular Dynamics simulation. 
We show that the theory of continuum elasticity based on nonlinear Fopple-Hencky equations is 
applicable on a circular suspended graphene flake. Our simulation results are plenty compatible to 
those results calculated by the nonlinear elasticity theory. We find the force-displacement curves 
in good agreement with the recent experimental measurements and conclude they are temperature 
independent. In addition, we find the vibration frequency for such a system and monitor that it 
behaves as a sinusoidal manner at small circular graphene size. 
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Introduction — Understanding of the mechanical proper- 
ties of graphene, a newly realized two-dimensional elec- 
tron system ^ is crucial in guiding its growth and 
applications. Booth et al Q have measured mechani- 
cal properties in a series of experiments on a suspended 
graphene and showed that graphene has an extraordi- 
nary stiffness which can supports an additional weight 
of many crystalline copper nano-particles. A nonlinear 
behavior of stress-strain response of free standing mono- 
layer graphene has been found by Lee et al [3] in pushing 
down a rigid tip of atomic force microscopy (AFM) which 
is placed over graphene. They determined that graphene 
is a very strong material as much strong as diamond. Fur- 
thermore, based on experimental measurements, a linear 
slope of force-displacement curves (FDCs) has been mea- 
sured. [1,0] 

The physical scales available to computer simulations 
strongly depend on the model. In general, the level of de- 
tail that the model encompasses reduces both the number 
of atoms and length of times that can be simulated. A 
sophisticate approach is the density functional theory 
where the electronic structure of the atoms is calculated 
self-consistently and it can generally treat only on the 
system incorporate in order of several hundred atoms. 
Classical Molecular Dynamics simulation (MD), on the 
other hand is a common method to do several physi- 
cal calculations for nano-indentation on various materi- 
als [1, Q to analysis systems containing huge number of 
particles with time scales on the order of nanoseconds. 
Accordingly, we use MD to study nano-indentation of 
a circular graphene flake (CGF) and to the best of our 
knowledge there is no any systematic study of simulated 
nano-indentation of the CGF. 

In experiment, since the hardness of AFM's tip is re- 
quired the diamond tip material is typically used. In 
general, geometrical properties of the three dimensional 
indenter on a bulk material are important either in ex- 
periments or simulations of nano-indentation. Q In order 
to check whether or not the existence of a geometric de- 
pendency, we examined two different structures, the di- 
amond structure with carbon atoms and a face-centered 




FIG. 1: (color online). Left panel: A snapshot of an in- 
denter is placed over a roughed graphene circular flake with 
cramped boundary conditions at room temperature. Right 
panel: A snapshot of circular graphene flake under influence 
of boundary tension before loaded by the indenter. 



cubic structure incorporating Platinum atoms, for an in- 
dentation simulation and found no indication of geomet- 
rical dependence in two dimensional graphene flake. In 
other words, our results indicate that the results are in- 
dependent of the geometry of indenter at least in the 
two structures considered for the indentation problem in 
graphene. 

The deformation of CGF under a central point load 
by carrying out molecular mechanics simulations has 
been investigated. In the recent work, authors used 
properly selected parameters and showed that the von 
Karman plate theory can provide a remarkable accurate 
prediction of the graphene sheet behavior under linear 
and nonlinear bending and stretching. In the present 
letter, we have used the theory of elasticity based on 
Fopple-Hencky equations for a CGF and shown that the 
theory is applicable to graphene sheets by comprising 
the results of theory with those results calculated by MD 
simulations. 

Theory and Model — We have used MD method to sim- 
ulate nano-indentation of a suspended CGF. Number 
of carbon atoms in our simulation varies from 18740 
to 110096 correspond to graphene surfaces with radius 
R = 12toi? = 30.1 nm, respectively. A rigidly 
clamped boundary condition has been employed. In 
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our simulations, the indenter is composed of 371 atoms. 
We have simulated the system at room temperature, 
(300i^) by employing Nos'e-Hoover thermostat. The 



Brenner's bond-order potential pjj, 112[ has been used 
for carbon-carbon interactions and Van der-Waals po- 
tential for the indenter-graphene interactions [13]. For a 
two-component system, as studied here, the parameters 
for the mixed interaction between two type of atoms can 
be estimated by the simple average suggested by Steel 
et al 



In the beginning of simulation, the lowest 
position of tip's atoms is considered a few angstroms, 
i.e. ^ 2.7 A above graphene sheet. In the load process, 
the indenter is pushed down slowly as 0.1 A over 5000 At 
which is equivalent to a velocity equal to 4 m/s. To avoid 
unphysical effects due to step time, the length of mov- 
ing indenter steps was chosen to be small with respect 
to the force cutoff length for the interatomic potential. 
The indenter atoms interact only with graphene atoms 
via the Lennrad- Jones interactions. The shape of the in- 
denter was chosen as a square-based pyramid as shown 
in Fig. 1. 

We have also applied the theory of elasticity f]3\ for 
the CGF under concentrated force at large deflection, z. 
Note that a typical size of system is much larger than the 
deflection values, and in addition the thickness of flake is 
much smaller than deflection amounts. 

The energetics considerations of a graphene in the limit 
of large deflections include both bending and stretch- 
ing energies, p^] The condition of minimum energy for 
graphene flake yields the Foppl equation. The solutions 
of equations can be given by using the Hencky transfor- 
mations. The governing equations in planar-polar coor- 
dinate are written as follow 
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where r is the radial position, R is the radius of CGF, 
z{r) is the deflection, t is the thickness of graphene ap- 
proximated by 1 A, is the Young's modulus, is the 
radial stress of flake and F is the concentrated load on 
the flake. We are interested on rigidly clamped bound- 
ary condition or fixed boundary condition in the absent 
of residual stress, i.e^ z = at r = and it impose the 
following equation [l5[ 
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where u is Poisson's ratio. Expressions given by Eq. ([T|) 
are nonlinear equations. At large deflection value, the 
set of equations can be solved analytically [16] with 
considering Eq. ([2j) and the concentrated force leads 
F = f^Q^C^ where G{i') is a complicated function 

of the Poisson ratio. Here C = I C I where ( = —z{0)k 
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FIG. 2: (color online). Top: Load force as a function of 
displacement. The fitting function given by Eq. (|3]) for three 
different sizes are shown by solid curves. Filled symbols are 
experimental data, Ref In the Inset: Load and unload 
( the indenter pulled up from maximum deflection position 
at C = 27 nm ) force as a function of displacement for the 
graphene flake with R = 12 nm. Bottom: Distribution of 
stress on graphene flake in arbitrary units when the indentor 
was at C = 27 nm. 



is graphene deflection at the center of CGF, r = in 
the k direction. G{iy) has almost a linear behavior with 
respect to Poisson's ratio of desired systems and mostly 
varies from 0.9 to 0.65 for given v ^ 0.0 to 0.5. Surprising 
enough, our simulation results confirm the concentrated 
force behavior which we will discuss that in discussion 
Section. 

The mechanical properties of the free standing 
graphene can be probed by the tip of AFM. The ex- 
perimental results have obtained for deflections greater 
than lOnm. Recently, the FDCs have been measured by 
Lee et al 0] and they have shown that the FDCs can be 
approximated by a simple polynomial function having a 
linear term due to tension effect and a cubic term due to 
bending effect. 



F = aC^bC, 



(3) 



and thus the bending term will dominate at large deflec- 
tion. The radius dependence of Hooke coefficient can be 
determined from the theory of elasticity 0] in which the 
Euler-Lagrange (EL) equation of minimization of the free 
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FIG. 3: (color online). Values of parameter 6 as a function 
of graphene circular size denoted by symbols. The solid line 
shows a 1/B? function. 



energy is —(jV'^z{~x) + i<iV^z{~x) = F5{~x — xq) where 
the (J and n are respectively the tension and the bending 
rigidity and = ix^y). Accordingly, the Hooke coeffi- 
cient of Eq. ([3]) has a complicated dependence on R. We 
have z{q) = Fjiaq^ + Kq^^ from the EL equation and it 
leads that a ~ — 47rcr/ ln(cr/ {gB? + /^)) . In this work, we 
therefore try to find a relationship between the parameter 
h with the radius of CGF. 

Discussions — Figure 1 (left panel) shows a snapshot of 
simulated CGF with an indenter placed over graphene 
(the arm of tip is shown only for clarity). The z- 
component of forces from graphene atoms on indenter are 
calculated by summing over total reaction forces. Note 
that the deflection or displacement of CGF is equiv- 
alent to the indenter displacement in each step. It is 
essential to compare the simulation results of FDCs with 
those measured in experiments for the cramped CGF. 

Figure 2 (top) shows the variation of applied load 
at r = as a function of graphene's deflection in z- 
direction. Corresponding values for parameters a and b 
have been obtained by fitting the numerical results with 
the expression given by Eq. ([3]). The size dependence 
of FDCs is determined in the figure. It is worthwhile 
to note that the sample size in the experiment possess- 
ing radius 500 nm is much larger than the one in our 
simulation samples. The FDC results for R = 30.1 nm 
where a = 0.108 N/m and b = 2.03 x lO^^N/m^, are very 
close to the data of experiments. Note that there are 
defects and impurities in a real sample and those would 
reduce the FDCs values. Apparently, studying the in- 
dentation properties for either a tiny graphene sheet or a 
tiny tip should be complicated in experiment. Obviously, 
the physics of indentation for tiny graphene sheet can be 
understood based on our MD simulations. 

One important physics which can be extracted from 



our simulation is unload FDCs. In the unload case where 
graphene was already loaded and reached to a maximum 
deflection value where C = 27 nm and then the indenter is 
pulled up, we have found that the force result in this case 
is the same as it in loaded case. In the inset of Fig. 2, the 
force results as a function of displacement for both load 
and unload cases are shown. It is physically known that 
the load and unload FDCs would be different when a bulk 
material is used. The physical reason is the deformation 
of structure make in the unload process. The distribu- 
tion of stress [13] on the CGF which depends on a mass 
density. Young's modulus and Poisson's ratio is shown in 
Fig. 2 (bottom) as well. The large stresses generated near 
the indenter at the center and decays fast as a function 
of radial value since the interaction of carbon-carbon is 
nonharmonic in the model interaction. In addition, the 
distribution of stress is no longer uniform on the surface 
due to the thermal fluctuations. 

In order to get quantitative results, we compare the 
concentrated force calculated from the theory and sim- 
ulation results fitted by Eq. (|3]) to describe the radial 
dependence of parameter b. The parameter b is given by 
b = 4^py;^- We have obtained the value of b by sim- 
ulating many CGF sizes and results are shown in Fig 3 
as a function of R. The value of b decreases by increas- 
ing R and it is fitted very well to 1/R^ function ( the 
slop of line in log- log plot is quite close to -2). It ob- 
viously means that the theory of elasticity is applica- 
ble for graphene flakes. Furthermore, we obtained that 
7rEt/G{jy) ^ 4 X 165 (N/m) using the fitting coefficient. 
Using approximated values for t ^ 1.0 A and G ^ 0.8, 
the Young modules can be obtained and it is about 1.65 
Terra Pascal which is very close to the value predicte d by 
recent experiments and theoretical calculations 0,11, [l3|. 
Moreover, the a and are respectively obtained about 
0.16 N/m and 0.8 eV from the fitting values of Hooke 
coefficient which are very close to values have been pre- 
dicted for graphene. [1, 

Interestingly, the mechanical vibration in suspended 
few-layer graphene has been illustrated in experiment 
by using a scanning prob microscopy and found a res- 
onance frequency about 31 MHz. [20| From the simula- 
tion point of view, when the indenter is pulled upward 
from its initial position there will be a position, {z{R)k) 
where graphene sheet separates from the indenter (af- 
ter the range of interaction between the carbon atoms 
and the indenter atoms) and the surface vibrates with 
high frequency. We have obtained the size dependence 
of frequency for one layer graphene and u = 32.0 ± 0.1, 
19.0 ± 0.3 and 15.5 ± 0.45 in units of GHz for R = 12, 15 
and 18.4 nm, respectively shown that v decreases by in- 
creasing R. We monitor that the vibration frequency be- 
haves as a sinusoidal manner at small CGF size, (12 nm). 
Our finding regard to the frequency at small graphene 
size would be verified by experiments. 

In the other hand, graphene under influence of bound- 
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FIG. 4: (color online). Load force as a function of displace- 
ment for both CGF (red triangles) and GIBT (blue circles). 



ary tension (GIBT) has been studied in several experi- 
ments. [1,[22[ It would be worth to explore the FDCs be- 
havior of GIBT. To this purpose, we generated a GIBT 
by compressing the boundary of CGF about 1.2% of its 
initial diameter and then the indenter over the GIBT is 
loaded with cramping boundaries condition. 

We have shown a snapshot of a GIBT in Fig. 1 (right 
panel). The out-of plane height on the surface is around 
18 A which is much bigger than the CGF where extra 
tensions are zero. The results of FDCs for the both cases 
have been shown in Fig. 4 for = 18.4 nm. Apparently, 
the force value of GIBT is smaller than one calculated 
for the CGF. The reason is that the number of atoms in 
the board range of interaction potential in the GIBT is 
larger than the one in the CGF and hence the attractive 
part of total force increases. Accordingly, the total force 
in this case might be smaller than the one calculated in 
the CGF. Consequently, the effect of boundary tension 
is essential for studying FDCs in comparison with the 
experiment measurements. 

In summary, our simulation results show that the 
FDCs is almost independent of ordinary temperatures ( 
up to 400 K) and they are close to the experimental data. 
Moreover, we found that the force on indenter in the CGF 
is larger than the GIBT. Importantly, the vibration fre- 
quency of the CGF in order of GHz is calculated and 
shown that it decreases by increasing the CGF sizes and 
essentially a particular vibration having a resemblance to 
a sinusoidal behavior occurs at small R. Eventually we 
would like to emphasize that our numerical results show 
that the continuum elasticity theory can be used for the 
indentation problem in graphene. 
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